Analyses of swisscom data
Grid preparation
Ancillary data
Municipality boundaries
st_gg21 <- st_read("data-raw/swisstopo/swissboundaries3d_2021-07_2056_5728/SHAPEFILE_LV95_LN02/swissBOUNDARIES3D_1_3_TLM_HOHEITSGEBIET.shp",
as_tibble = TRUE) %>%
st_zm(drop = TRUE, what = "ZM") %>%
# move to LV03
st_transform(21781) %>%
# filter(! BFS_NUMMER %in% c(2391, 5391, 5394)) %>%
# exclude lakes
# filter(OBJEKTART != "Kantonsgebiet") %>%
# exclude FL & enclaves
filter(ICC == "CH") %>%
select(BFS_NUMMER, NAME, KANTONSNUM, GEM_TEIL) %>%
rename(GMDNR = BFS_NUMMER,
GMDNAME = NAME,
KTNR = KANTONSNUM) %>%
arrange(GMDNR)
write_rds(st_gg21, "data/swisstopo/st_gg21.Rds")STATPOP offset
Could be used to define offset for st_make_grid
statpop_20 <- read_delim("data-raw/BfS/ag-b-00.03-vz2020statpop/STATPOP2020.zip",
delim = ";", escape_double = FALSE, trim_ws = TRUE)[] %>%
mutate_all(as.integer)
# statpop_20 %>% slice(1:10) %>% View()
write_rds(statpop_20, "data/BfS/statpop_20.Rds")offset_bfs <- read_rds("data/BfS/statpop_20.Rds") %>%
summarise(min_x = min(X_KOORD),
min_y = min(Y_KOORD)) %>%
st_as_sf(coords = c("min_x", "min_y"),
crs = 21781,
remove = FALSE)swisscom offset
Using two communities that are furthest away in southerly and westerly directions:
Chansy to define x
- Postal code 1284
- SFOS number 6611
Chiasso to define y
- Postal code 6830
- SFOS number 5250
Tile definitions were pulled from API using query_postal_codes_heatmaps_api.py.
Points of grid were defined using lower left corner coordinates.
WGS84 coordinates were transformed to LV03.
chansy_grid <- fromJSON("data/swisscom/chansy_grid.json")
chansy_grid <- flatten(chansy_grid$tiles) %>%
as_tibble()
chansy_grid_sf <- chansy_grid %>%
st_as_sf(coords = c("ll.x", "ll.y"),
crs = 4326,
remove = FALSE) %>%
st_transform(21781) %>%
mutate(x = st_coordinates(.)[, 1],
y = st_coordinates(.)[, 2])chiasso_grid <- fromJSON("data/swisscom/chiasso_grid.json")
chiasso_grid <- flatten(chiasso_grid$tiles) %>%
as_tibble()
chiasso_grid_sf <- chiasso_grid %>%
st_as_sf(coords = c("ll.x", "ll.y"),
crs = 4326,
remove = FALSE) %>%
st_transform(21781) %>%
mutate(x = st_coordinates(.)[, 1],
y = st_coordinates(.)[, 2])Example of grid
Coordinates of lower left corner were then obtained by getting minimum values
offset_swisscom <-
bind_rows(
chansy_grid_sf %>%
st_drop_geometry() %>%
summarise(min_x = min(x),
min_y = min(y)) ,
chiasso_grid_sf %>%
st_drop_geometry() %>%
summarise(min_x = min(x),
min_y = min(y))
) %>%
summarise(min_x = min(min_x),
min_y = min(min_y)) %>%
# needs rounding - transfrom error?
mutate(min_y = round(min_y)) %>%
st_as_sf(coords = c("min_x", "min_y"),
crs = 21781,
remove = FALSE)Note small difference from BfS derived minimums (in black)!
Grid
Generate whole country
100m grid, with lower left corner defined using STATPOP cells.
Also, adding ID for each cell (based on row number).
country <- st_gg21 %>%
st_make_grid(cellsize = 100,
offset = c(offset_swisscom$min_x, offset_swisscom$min_y),
square = TRUE) %>%
st_sf() %>%
st_cast("POLYGON") %>%
mutate(ID1 = row_number()) %>%
relocate(ID1)
# st_write(country, "data/grid/country.gpkg")
write_rds(country, "data/grid/country.Rds")Generate cities
Clipped / selected using swisstopo municipality data.
Added second, city specific ID and area.
Function with clipping
# why st_intersection
# https://stackoverflow.com/questions/62442150/why-use-st-intersection-rather-than-st-intersects
city_grid_clip <- function(municipality){
city <- country %>%
st_intersection(st_gg21 %>%
filter(GMDNR == municipality)) %>%
mutate(ID2 = row_number(),
area = st_area(.)) %>%
relocate(ID1, ID2, area)
return(city)
}Function without clipping
# solution from
# https://stackoverflow.com/questions/57014381/how-to-filter-an-r-simple-features-collection-using-sf-methods-like-st-intersect
city_grid <- function(municipality){
city <- country %>%
filter(st_intersects(geometry,
st_gg21 %>% filter(GMDNR == municipality),
sparse = FALSE)) %>%
mutate(ID2 = row_number()) %>%
relocate(ID1, ID2)
return(city)
}# or even simpler, without dplyr
# could be generalized with arguments to provide full grid and GMDE shape
city_grid <- function(municipality){
city <- country[st_gg21 %>% filter(GMDNR == municipality), ] %>%
mutate(ID2 = row_number()) %>%
relocate(ID1, ID2)
return(city)
}Zuirich
# zurich_clip <- city_grid_clip(261)
zurich <- city_grid(261)Geneva
# geneva_clip <- city_grid_clip(6621)
geneva <- city_grid(6621)Basel
# basel_clip <- city_grid_clip(2701)
basel <- city_grid(2701)Lausanne
# lausanne_clip <- city_grid_clip(5586)
lausanne <- city_grid(5586)Bern
bern_clip <- city_grid_clip(351)
bern <- city_grid(351)Offset
Using two communities to get min x and y coordinates: Chansy & Chiasso
# chansy_clip <- city_grid_clip(6611)
# st_write(chansy_clip, "data/grid/chansy_clip.gpkg")
# write_rds(chansy_clip, "data/grid/chansy_clip.Rds")
chansy <- city_grid(6611)
# st_write(chansy, "data/grid/chansy.gpkg")
write_rds(chansy, "data/grid/chansy.Rds")
# chiasso_clip <- city_grid_clip(5250)
# st_write(chiasso_clip, "data/grid/chiasso_clip.gpkg")
# write_rds(chiasso_clip, "data/grid/chiasso_clip.Rds")
chiasso <- city_grid(5250)
# st_write(chiasso, "data/grid/chiasso.gpkg")
write_rds(chiasso, "data/grid/chiasso.Rds")Example
Municipality Bern coverage:
Clipped grid:
Vs. unclipped one:
Note, when using clipped grid, there are few examples of grid cells with very small area.
These could be excluded if no population was found there?
ID1 ID2 area GMDNR GMDNAME KTNR GEM_TEIL
1 4259841 401 0.1571798 [m^2] 351 Bern 2 0
2 4332927 2876 0.1766709 [m^2] 351 Bern 2 0
3 4423647 5188 0.2123485 [m^2] 351 Bern 2 0
4 4339904 3098 0.2957572 [m^2] 351 Bern 2 0
5 4273711 702 0.7733263 [m^2] 351 Bern 2 0
6 4245805 139 1.0971039 [m^2] 351 Bern 2 0
7 4245786 123 1.1456480 [m^2] 351 Bern 2 0
8 4420099 5082 1.9921709 [m^2] 351 Bern 2 0
9 4291217 1374 2.4771009 [m^2] 351 Bern 2 0
10 4273760 735 3.2587531 [m^2] 351 Bern 2 0
Area estimates
Bern - hourly
Number of tiles for area of Bern
(tiles <- nrow(bern))[1] 5487
Number of requests needed
(requests <- ceiling(tiles / 100))[1] 55
Price for complete grid
(price <- requests * 0.01)[1] 0.55
Price for 24h
(price_24h <- price * 24)[1] 13.2
Price for 4w hourly
(price_4w <- price_24h * 7 * 4)[1] 369.6
Price for year hourly
(price_year <- price_24h * 365)[1] 4818
Price for 4w daily
(price_4w_d <- price * 7 * 4)[1] 15.4
Price for year daily
(price_year_d <- price * 365)[1] 200.75
5 cities - hourly
Number of tiles for area of Bern
(tiles <- nrow(zurich) +
nrow(geneva) +
nrow(basel) +
nrow(lausanne) +
nrow(bern))[1] 24116
Number of requests needed
(requests <- ceiling(tiles / 100))[1] 242
Price for complete grid
(price <- requests * 0.01)[1] 2.42
Price for 24h
(price_24h <- price * 24)[1] 58.08
Price for 4w hourly
(price_4w <- price_24h * 7 * 4)[1] 1626.24
Price for year hourly
(price_year <- price_24h * 365)[1] 21199.2
Price for 4w daily
(price_4w_d <- price * 7 * 4)[1] 67.76
Price for year daily
(price_year_d <- price * 365)[1] 883.3